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Abstract The present work studies the isotropic and homogeneous turbulence for incompressible 
fluids through a specific Lyapunov analysis, assuming that the turbulence is due to the bifurcations 
associated to the Navier-Stokes equations. 

The analysis consists in the calculation of the velocity fluctuation through the Lyapunov analysis 
of the local deformation and the Navier-Stokes equations and in the study of the mechanism of the 
energy cascade from large to small scales through the finite scale Lyapunov analysis of the relative 
^ ■ motion between two particles. 

The analysis provides an explanation for the mechanism of the energy cascade, leads to the closure 
^ . of the von Karman-Howarth equation, and describes the statistics of the velocity difference. 

O ' Several tests and numerical results are presented. 
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1 Introduction 
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t^hj- , In this work a novel procedure based on a specific Lyapunov analysis is presented for studying the 

incompressible isotropic and homogeneous turbulence in an infinite domain. The analysis is mainly 
motivated by the fact that in turbulence the kinematics of the fluid deformation is subjected to bifur- 
cations pQ and exhibits a chaotic behavior and huge mixing [2], [3], resulting to be much more rapid 
than the fluid state variables. This characteristics implies that the accepted kinematical hypotheses for 
deriving the Navier-Stokes equations could require the consideration of very small length scales and 
times for describing the fluid motion [3] and therefore a very large number of degrees of freedom. 
As well known, other peculiar characteristics of the turbulence are the mechanism of the kinetic energy 
^ ■ cascade, directly related to the relative motion of a pair of fluid particles [5], [B], [7], [5] and respon- 

sible for the shape of the developed energy spectrum, and the non-gaussian statistics of the velocity 
difference. 

This energy spectrum can be calculated through a proper closure of the von Karman-Howarth 
equation (see Appendix) or of its Fourier Transform [7] , [5] . The equation describes the evolution of 
the correlation function / of the longitudinal velocity u r , and depends upon the term K (see Appendix), 

Department of Mechanics and Aeronautics 
University "La Sapienza", Rome, Italy 
via Eudossiana, 18 
Tel.: +39-06-44585268 
Fax: +39-06-4881759 
E-mail: dedivitiis@dma.dma.uniromal.it 



2 



Nicola de Divitiis 



directly related to the longitudinal triple-velocity correlation k. This latter, due to the inertia forces, 
does not change the kinetic energy of the fluid and satisfies the detailed conservation of energy [H] 
which states that the exchange of energy between wave-numbers is only related to the amplitudes of 
these wave-numbers and of their difference [5] . 

Various authors (see for instance [10], [11], [12]), propose, for k, the following diffusion approxima- 
tion 

k = ^f in 

u or 

where r and D = D(r) are the separation distance and the turbulent diffusivity, whereas u 2 = (uiUi}/3 
represents the longitudinal velocity standard deviation. This implies that the closed von Karman- 
Howarth equation is a parabolic equation in any case (also for v =0). 

To the author knowledge, Hasselmann in 1958 [IU] was the first that proposed a link between k 
and /, using a simple model which expresses k in function of the momentum convected through the 
surface of a spherical volume. His model incorporates a free parameter and expresses D(r) by means 
of a complex expression. 

Another closure model in the framework of Eq. ([1} was developed by Millionshtchikov [TT]. There, 
the author assumes that D(r) = k\u r, where k\ is an empirical constant. Although both the models 
describe two possible mechanisms of the energy cascade, in general, do not satisfy some physical 
conditions. For instance, the model of Hasselmann does not verify the continuity equation for all the 
initial conditions, whereas the Millionshtchikov's model gives, according to Eq. (jTJ) , an absolute value 
of the skewness of du r /dr, in contrast with the several experiments and with the energy cascade [8]. 

More recently, Oberlack and Peters [T2] suggested a closure model where D is in terms of /, i.e. 
D{r) = A^r — f and is a constant parameter. The authors show that this closure reproduces the 
energy cascade and, for a proper choice of k2, provides results [T2] in agreement with the experimental 
data of the literature. 



In general, Eq. (TJJ represents models of diffusion approximation based on the assumption that the 
turbulence can be represented by an opportune diffusivity which varies with r [5] . As a consequence of 
Eq. (TJJ, K contains a term proportional to d 2 f/dr 2 which can occur only if the inertia forces include 
stochastic external terms, independent from the fluid state variables [13] and not present in the classical 
formulation [7J, [H]- For this reason the models based on Eq. (TT]) are really phenomenological closure 
ofEq. ([7H]). 

Although several other works on the von Karman-Howarth equation were written [T3] , [TS] , [TB] , [T7J , 
to the author's knowledge, a theoretical analysis based on basic principles which provides a physical- 
mathematical closure of the von Karman-Howarth equation and the statistics of Au r has not received 
due attention. Therefore, the objective of the present work is to develop a theoretical analysis based 
on reasonable physical conjectures which allows the closure of the von Karman-Howarth equation and 
the determination of the statistics of Au r . 

Of course, besides the von Karman-Howarth equation, there are some other approaches, such as the 
direct numerical simulation based on solving the Navier-Stokes equations and the experiments, which 
are not considered in the present work. 

The present work only considers the possibility to obtain the fully developed homogeneous-isotropic 
turbulence in a given condition and does not analyze the intermediate stages of the turbulence. The 
study assumes that the fluctuations of the fluid state variables are the result of the bifurcations of the 
Navier-Stokes equations. In section[5] we present a qualitative scenario of these bifurcations which leads 
to the onset of the turbulence. These bifurcations, defined by means of the fixed points of the velocity 
field and the Navier-Stokes equations, allows a rough estimation of the critical Reynolds number based 
on the Taylor scale. After, in section [3] the velocity fluctuation is studied through the kinematics of 
the local deformation and the momentum equations. These latter are expressed with respect to the 
referential coordinates which coincide with the material coordinates for a given fluid configuration [¥] , 
whereas the kinematics of the local deformation is analyzed with the Lyapunov theory. The choice 
of the referential coordinates allows the velocity fluctuations to be analytically expressed in terms 
of the Lyapunov exponent of the local fluid deformation. The section 2] deals with the study of the 
velocity difference between two fixed points of the space. This is analyzed with an opportune finite 
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scale Lyapunov theory studying the motion of the particles crossing the two points, in the finite scale 
Lyapunov basis. This basis is formally obtained through the orthonormalization method of Gram- 
Schmidt of the finite scale Lyapunov vectors. These latter, also known as Bred vectors, were first 
introduced by Toth and Kalnay [18 to study the evolution in the time of a nonlinear perturbed model 
subjected to an initial finite perturbation. The choice of such vectors, whose properties are related to 
the classical Lyapunov vectors , is revealed to be an usefull tool for representing the relative motion 
of two particles crossing two given points of the space. 

The present analysis postulates that the motion of such Lyapunov basis and that of the fluid with 
respect to the same basis, are completely statistically uncorrelated. This crucial assumption arises 
from the condition of fully developed turbulence. The study leads to the closure of the von Karman- 
Howarth equation 7 and gives an explanation of the mechanism of the kinetic energy transfer between 
length scales. 

The obtained expression of K is the result of this assumption and does not correspond to a diffusive 
approximation with model free parameters. Its mathematical structure is in terms of / and df / dr and 
satisfies the conservation law which states that the inertia forces only transfer the kinetic energy [7] , [5] . 
This expression of K corresponds to a first order term which makes the closed von Karman-Howarth 
equation, a nonlinear partial differential equation of the first order in r when v = 0. The main asset of 
the proposed closure on the other models is that it has been derived from a specific Lyapunov theory, 
with reasonable basic assumptions about the statistics of the velocity difference. 

Furthermore, the statistics of the velocity difference is studied in section [7] with the Fourier analysis 
of the velocity fluctuations, and an analytical expression for the velocity difference and for its PDF is 
obtained in case of isotropic turbulence. This expression incorporates an unknown function, related to 
the skewness, which is identified through the obtained expression of K. This velocity difference also 
requires the knowledge of the critical Reynolds number whose estimation is made in the section [3] 

Finally, the several results obtained with this analysis are compared with the data existing in the 
literature, indicating that the proposed analysis adequately describes the various properties of the fully 
developed turbulence. 



2 Bifurcations 

This section qualitatively describes the route toward the turbulence by means of the bifurcations of 
the Navier-Stokes equations and provides an estimation of the critical Taylor scale Reynolds number, 
assuming that the turbulence is fully developed, homogeneous and isotropic. 

The velocity field u = u(x, t) of a viscous and incompressible fluid, measured in the reference frame 
5i, satisfies the Navier-Stokes equations 

V* • u* = 

du* ( 1 , \ ( 2 ) 

= - u*V*u* + vv V* V 

dt* \ F Re 

Into Eq. ©, u* = u/Z7, t* = tU/L, x* = x/L, p* — p/pU 2 , and Re = UL/i>, where U and L are 
assigned velocity and length, respectively. The pressure p can be eliminated by taking the divergence 
of the momentum equation [§]. The velocity field, starting from the unique initial condition u(x, 0) 
which does not depend on the Reynolds number, will depend upon Re by means of its time evolution 

^=F(xV*;i?e) (3) 

where F(x*, t*; Re) represents the right-hand-side of the momentum Navier-Stokes equations calculated 
for u* = u*(x*,f*). 

Consider now the velocity field at t = 0, and the fixed points X of Eq. © which, by definition, 
satisfy du* /dt* = [50]. Increasing the Reynolds number, X will vary according to Eq. ©, which can 
be expressed through the implicit function theorem [20] 

r^ e <9F 

X = X - / VF- 1 — dRe (4) 

jReo dRe 
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where Re plays the role of the control parameter and Xo is the fixed point calculated at Re — ReQ, 
for t = 0. The location of these points will depend on the momentum Navier-Stokes equations and on 
the mathematical structure u(x, 0) [50] ■ Therefore, Re influences the distribution of X in the space. 

According to the literature [21], [35], [53], [23] and to the characteristics of the diverse kinds of 
bifurcations, we assume the following qualitative scenario: 

For small Re, the viscosity forces are stronger than the inertia ones and make F an almost smooth 
function of X. When the Reynolds number increases, as long as the Jacobian VF is nonsingular, X 
exhibits smooth variations with Re, whereas at a certain Re, this Jacobian becomes singular (det (VF) 
= 0). This can correspond to the first bifurcation, where at least one of the eigenvalues of VF crosses 
the imaginary axis and X appears to be discontinuous with respect to Re [50] . 

Figure [JJ shows a scheme of bifurcations at t — 0, where the component X of X is reported in terms 
of Re. Starting from Reo, the diagram is regular, until Rep, where the first bifurcation determines 
two branches whose maximum distance is AX p. AX and ARe give, respectively, a length scale of the 
velocity field at the current value of Re, and the distance between two successive bifurcations. After P, 
Eqs. fl3J and (J3J) do not indicate which of the two branches the system will choose, thus a bifurcation 
causes a lost of informations with respect to the initial data [25) . That is, very small variations on the 
initial condition or very little perturbations, are of paramount importance for the choice of the branch 
that the system will follow [25 . This is the situation of the bifurcations at t = 0. 

As the time increases, the position of the fixed points vary and so also the bifurcations, and far 
from the initial condition, one observes a developed motion, where the bifurcations can continuously 
vary with the time. Therefore, the bifurcations map changes with the time, and, if Re is high enough, 
the length scales are continuously distributed. In this case the energy spectrum can be continuous and, 
according to the theory [215], the velocity behaves as a chaotic function of t and x. There, the scales of 
Taylor and Kolmogorov are supposed to be assigned steady quantities. 

To describe the road to turbulence, observe that the average distances l n = (AX n ) depend on Re 
through Eqs. (J3J, and are here approximated by 20 

l n = 4h (5) 
a" 1 

where a » 2, [2TJ and the average is calculated on the time. Equation ([5]) is supposed to describe 
the route toward the chaos and is assumed to be valid until the onset of the turbulence. There, the 
minimum for L can not be less than the Kolmogorov scale I = (v 3 /e) 1/A [JJ, [2S] where h gives a good 
estimation of the correlation length of the phenomenon [20], [25] which, in this case is the Taylor scale 
At- Thus, i < l n < At, and 

£=^h (6) 




Fig. 1 Map of the bifurcations at an assigned time. 
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where N is the number of bifurcations at the beginning of the turbulence. Equation ([6]) gives the 
connection between the critical Reynolds number and N. In fact, the characteristic Reynolds num- 
bers associated to the scales I and At are Rk — Iuk/v = 1 and R\ = Xtu/v, respectively, where 
ur = (^e) 1 ^ 4 is characteristic velocity at the Kolmogorov scale, and u = y/ (ujUj) /3 [8 . For isotropic 
turbulence, these scales are linked each other by [5] 



\ T /t=15 1 '*y/R* (7) 
In view of Eq. ([6]), this ratio can be also expressed through N, i.e. 

a N-l = 15 l/4^ (8) 



Assuming that a is equal to the Feigenbaum constant (2.502...), the value R\ ~ 1.6 obtained for N = 
2 is not compatible with At which is the correlation scale, while the result R\ ~ 10.12, calculated 
for N = 3, is an acceptable minimum value for R\. The order of magnitude of these values can be 
considered in agreement with the various scenarios describing the roads to the turbulence [35], [2Tj . 
[23] . [2~i] . and with the diverse experiments [37], [25], [5S] which state that the turbulence begins for 
N > 3. Of course, this minimum value for R\ is the result of the assumptions a ~ 2.502, l\ ~ At, 
I n — £ and of approximation ([5]) . 



3 Lyapunov analysis of the velocity fluctuations 

In this section, the velocity fluctuations caused by the bifurcations of Eqs. ([5]) [T], are studied through 
the Lyapunov analysis of the kinematic of the fluid strain, using the Navier-Stokes equations. 

Starting from the momentum Navier-Stokes equations written in a frame of reference 3? 
du k du k , 1 dT kh 

—zrr — —t, — u h H q — (9 

at axh p oxh 

consider the map x : xo — > x, which is the function that determines the current position x of a fluid 
particle located at the referential position Xrj [1] at t = to- Equation Q can be written in terms of the 
referential position Xo [I] 
du k ( du k t 1 dT kh \ dx 0p 



Uh + 1^±1L (10 ) 
at \ dxQ P p dxQ P J dx h 

where Tkh represents the stress tensor. The Lyapunov analysis of the fluid strain provides the expression 
of this deformation in terms of the maximal Lyapunov exponent 

*L w e A(t-t ) 

9x0 

where A = max(yli, /I2, A3) is the maximal Lyapunov exponent and A4, (i = 1, 2, 3) are the Lyapunov 
exponents. Due to the incompressibility, Ai + A% + A3 = 0, thus, A > 0. 

The momentum equations written using the referential coordinates allow the factorization of the 
velocity fluctuation and to express it in Lyapunov exponential form of the local fluid deformation. If 
we assume that this deformation is much more rapid than dTkh/dxo P and duk/dxo p Uh, the velocity 
fluctuation can be obtained from Eq. (Tit)]) , where dTkh/ dx$ p and duk/dxopUh are supposed to be 
constant with respect to the time 

Uk ~ a v dx 0p Uh pdx 0p ) t=to ~ a\ at y t=iD 

This assumption is justified by the fact that, according to Truesdell @], 8T k h / 9xq p — duk/dxo p Uh is a 
smooth function of t -at least during the period of a fluctuation- whereas the fluid deformation varies 
very rapidly in the vicinity of a bifurcation according to Eq. (jlip . This implies that, in proximity of 
bifurcations, the Lyapunov basis of orthonormal vectors Ea = (e'^e^e^) [30j associated to the strain 
(fTTjl rotates very quickly with respect to with an angular velocity oja, whereas the modulus of the 
fluid velocity fluctuation, measured in Ea increases with a rate ps e^* - ' ). Since A is related to the 
maximal eigenvalue of (Vu + Vu T )/2, according to the fluid kinematics [2], [3], [3T] \u>a\ = 0(A). 
Therefore, the fluid velocity fluctuation measured in 3?, varies very quickly because of the combined 
effect of the exponential growth rate and of the rotations of Ea with respect to 3?. 

Note that, since 9x/9xo is supposed to vary much more fastly than the fluid state variables, as long 
as t — to does not exceed very much the Lyapunov time 1 /A, the Lyapunov vectors almost coincide 
with the eigenvectors of Vu and A is nearly the maximum eigenvalue of (Vu + Vu T )/2 |20j . 
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Fig. 2 Scheme of the relative motion of two fluid particles 



4 Lyapunov analysis of the relative motion 

In order to investigate the mechanism of the energy cascade, the relative motion between two fluid 
particles is here studied with the Lyapunov analysis. 

Consider now two fixed points of the space, X and X' (see Fig. [2]) with r = X' — X, where r = |r| is 
the separation distance, and the motion of the two fluid particles which at the time to, cross through 
X and X'. The equations of motion of these particles are 

§ =u(M), ^ =»(*',«) w 

where u(x, t) and u(x',t) vary with the time according to the Navier- Stokes equations, f = x' — x is 
the relative position between the two particles, therefore r(to) = r 

The Lyapunov analysis of Eqs. (| 13[) leads to the determination of the maximal finite scale Lyapunov 
exponent A of Eqs. © and ([TSf 



A(r) = lim — / — f dt (14) 

where A(0) = A, and r represents the scale associated to A. To define the finite scale Lyapunov basis 
of Eqs. (fT5|) . consider the system 

f=„(M) (15) 

^ = u(x + r,t)-u(x,i) (16) 
dt 

The finite scale Lyapunov vectors ri, Y2 and r3 are defined as the solutions of Eq. (flo]) , starting 
from a given initial condition r^(0), (i =1, 2, 3) jTS], where x and u vary according to Eqs. (IT5|) and 
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([2]), respectively. Such vectors, also known as Bred vectors, are, by definition, related to the classical 
Lyapunov vectors in such a way that, ri, r 2 and r3 tend to the classical Lyapunov vectors when 
r-j(O) — > 0, (i =1, 2, 3) [T5], [35]. The finite scale Lyapunov basis E\ = (ei,e 2 ,e 3 ) of Eqs. (fTB")) is then 
obtained by orthogonalizing ri, r 2 and r3 with the Gram-Schmidt method at each time [3U], |32j . 
This basis rotates with respect to 3? with an angular velocity — u>\, where, because of isotropy 
f=s A(r). The fluid velocity difference Av = (v[ — v%, v' 2 — V2, v' 3 — V3), measured in Ex, is 
expressed by the Lyapunov theory as 

v't — vi = Xi ?i, 1 = 1,2,3 (17) 

Into Eq. (|I7p . A; is the Lyapunov exponent associated to the direction fi, where 77, vi and v[ are, 
respectively, the components of x' — x and of u(X,i) and u(X',£) measured in E\. Av can be also 
written in terms of the maximal exponent A 

Av(r) = A(r)r + £ (18) 

where due to the other two exponents, is negligible with respect to Ar and makes Av(r) a solenoidal 
field, resulting |£| << |Ar| during the fluctuation. When C/Ar — > 0, f maintains unchanged orientation 
with respect to E\. The velocity difference Au, measured in 3? and expressed in Ex, is calculated when 
|r| = r 

zlu(f) = A(r)f + C + w A xf (19) 
This equation states that u)\ determines the lateral components of the velocity difference in E\. 

The longitudinal component of the velocity difference in -K, Au r is then 

Au r {v) = £ ■ QZ\u(r) (20) 

where Q = ((e,,)) is the rotation matrix transformation from E\ to 5i, eij are the components of in 
3fi, and ^ = (X' — X)/|X' — X| is the unit vector along the longitudinal direction. Observe that, A(r) 
and U3\ arise both from the Navier-Stokes equations and that, in line with Eq. (|T^|) . A(r) can change 
but its variations are much slower than u>\ [20]. Moreover, the Lyapunov vectors sweep a subdomain 
of their space of state which coincides with the physical space R 3 . The simultaneous variations of these 
quantities produce the velocity difference fluctuation observed in 3?, according to Eq. (l20l) . 



Although u>\ and Av arise from Eqs. (|15p and (|16D . these are not directly related, therefore, we 
suppose here that lj\ is statistically orthogonal to both v and v', with (uj) = 0. This is the crucial 
assumption of the present work which is justified by the condition of fully developed turbulence. This 
provides that 

(uxvpv' q n ) = 0, m, n > 0, p, q = 1, 2, 3 (21) 

where the angular brackets denote the average calculated on the statistical ensemble of all the pairs 
of particles which cross through X and X'. Again thanks to the fully developed turbulence, we also 
suppose that Q and u)\ are statistically orthogonal each other. This implies that 

« pe ^) = « P )(er g ), m,n>0, i,p,q = 1,2,3 (22) 

whereas the isotropy gives 

(u\hU\k) = g^A ' U\)Shk, (ehpCkq) = T^ShkSpq (23) 

where u)\ = ^2,^3), and e l = (en, e l2 , e i3 ), (i = 1,2,3), are expressed in 3?, and 

(uli) =a,A 2 (r), i = l,2,3 (24) 

Because of homogeneity, et^ = O(l) do not depend on r. 

Under these hypotheses, we now show the following equations 
(u r u' r 0Jxk) = (25) 

(u n u' n uJxk) = (u b u' b ujxk) = C k u 2 \(r)g(r) (26) 
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where u r , u n and u& are the components of u along f and along the two orthogonal directions n 
and 6, and C k is a proper constant of the order of unity. The function g(r) = (u n u' n )/u 2 is the lateral 
velocity correlation function (see Appendix) measured in E\ , that due to the homogeneity and isotropy, 
coincides with the lateral correlation function measured in 5i. 

Equation ([25]) is the direct consequence of Eq. ([2"T]) . In fact (u r u' r u>x) = (v r v' r uj\) — 0. 

To demonstrate Eq. (|26p . observe that u n and oj\ can be decomposed into n independent, identically 
distributed, stochastic variables which satisfy [33J 

(&>=0, (&&>=* h k (Zhtel)=™ hkj q (27) 

where ro/t/y = 1 if h = A; = j, else t*7/ifej = and q ^ 0. w>j is expressed as the sum of where, 
without loss of generality, the coefficients of the combination are assumed constant with respect to X 
and equal each other [33J . 

1 ™ 

WAi = A,A(r)~£)& (28) 
fe=i 

with A\ = Ai = A-^ = 0(1) constant parameters. The component m„ (and Ub) is also expressed as the 
linear combination of , whose coefficients are now functions of X as the consequence of the previous 
assumption, i.e. 

n 

u n (X)=uJ2Fk(X)Hk (29) 
fe=i 

Hence, the correlation function g(r) is in terms of F k and is determined through Eq. (|2"9"]l and (|27D 
[33] . putting X = 0. 

n 

<?(r) = ^F fc (0)F fe (r) ( 3 °) 
fe=l 

Therefore, (u n u' n uj\k) is calculated taking into account Eq. ((28]) . f|29[) . (l30t and ([27]), and, as a result, 
Eq. ([21)]) is achieved. 



5 Closure of the von Karman-Howarth equation 

Now, we present the closure of the von Karman-Howarth equation based on the analysis seen at the 
previous section. 

The term representing the inertia forces in the von Karman-Howarth equation satisfies the identity 
(l80l) (sec Appendix), which is here written as [7], [S] 



d d 

0^- (r k K) = — (uiUiiuk - u k )) (31) 

where, due to the isotropy if is a function of r alone [7]. The divergence of rK gives the mechanism 
of the energy cascade which does not depend upon the frame of reference [7] [5] . In order to determine 
the expression of K, Eq. ([3T]) is here written in E\. In view of Eq. ([T9"]) and taking into account that 
Uiii! i is also frame independent, one obtains the following equation 

d d 

— (f k K) = -— ((uiu' 2 X) f k + {(u.u'iLJx) x r) k + (uiu'iCk)) (32) 
or k dr k 

Since A is calculated with Eq. (Tl4l) . this is constant with respect to the statistics of and u[, thus 
(Xuiu'j) = A(witt-), and K is expressed as the general integral of Eq. ([32"]) 

Ki = — A (uiii't) f — {uiu^x) x f + s (33) 

Into Eq. (|33[) . s is the sum of a term due to £ plus an arbitrary solenoidal field arising from the 
integration of Eq. ([32l [53] . According to this analysis of Lyapunov, s is proportional to u 2 Xr and can 
be written in the form 

s = u 2 \(r)r s (34) 
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Substituting Eq. ([76]) (see Appendix) and Eq. ([34]) into Eq. ([33]), Kr is 

Kr = (Xu 2 {g - /) - 3Xu 2 g) r - (u^ux) X f + u 2 \(r)r s (35) 

where f(r) is the longitudinal velocity correlation function. 

Equation (1331) is made by three addends. In the first one of these, the part into the brackets is an 
even function of r which goes to zero as r — > oo and assumes the value — 3m 2 A(0) for r = 0. The second 
term, orthogonal to the first one, vanishes at r — and tends to zero when r — > oo. In this latter Uiu[ 
is expressed in E\ for sake of convenience 

(mulux) = (u r u' r u>x) + (u n u' n u)x) + (u b u' b u)\) (36) 

The first term at the RHS of Eq. ([33]) vanishes because of Eq. (|2~5]) , whereas the other ones are expressed 
by means of Eq. ([26]) . Therefore 

{mu'jUx) x f = 2u 2 X(r)g(r)c x f (37) 

where c = (d, C 2 , C a ) and C k = 0(1) are from Eq. (26]). 

In order to satisfy Eq. ([3S"|). the expression of So must be of the kind So = h(r)t + p(r)n, where, 
because of homogeneity and according to Refs. [8] and [35], h(r) and p(r) are even functions of r. 
Moreover, due to the isotropy and without lack of generality, h(r) — p(r) [5], [35], so Sq is 

s = /i(r)(t + n) (38) 

where t = f/r, n = (cx f)/|cx f|. 

To determine K and h(r), Eq. (1331 is projected along the directions f and n 

K = \u 2 (g-f)-3u 2 \g+^ (39) 

s • n = {uiU^oJx x r) • n (40) 
From Eq. gD]) 

h(r) = Hg(r) (41) 

where H = 2(c x f) ■ n/r = 0(1). As £/r — > 0, f maintains unchanged orientation in Ex, thus H is an 
invariant which has to be identified. The function K is determined with Eq. (|39[) 

K = Xu 2 (g-.f)+u 2 X(r)g(r)(H-3) (42) 

K(r) satisfies the conditions dK(0)/dr — and K(0) = jS], which represent, respectively, the 
homogeneity of the flow and the condition that the inertia forces do not modify the fluid kinetic 
energy. Since g(0) — f(0) — 1, this immediately identifies H = 3 and 

K = \ u 2 (g — f) (43) 

Due to the fluid incompressibility, / and g are related each other through g = f + l/2df /dr r (see Eq. 
d78l) . Appendix), leading to the expression 



K = -u 2 ^- X(r)r (44) 
2 dr w v ; 

This expression of K has been obtained studying the properties of the velocity difference in Ex ■ 

Equation ([4~4"]l states that, the fluid incompressibility, expressed by g — f ^ 0, represents a sufficient 
condition to state that K ^ 0. This latter is determined as soon as A is known. To calculate A, it is 
convenient to express An = u(x',t) — u(x, t) in 5ft, with |f| = r. Thus, Au r is first expressed in terms 
of f and Av through Eq. ([2"0"]) (Au r = £ ■ QZ\u), then its standard deviation is calculated assuming 
that C — 0, (i.e. Z\u = Af + u>x x f) and taking into account that Q and ljx satisfy Eq. ([2"3]) , with (ojx) 
= 0: 

{(Au r ) 2 ) = J2i jkiJ^pq r s (6£p(A 2 eije pg )fjfg + S,iS, P {Xe lj e pq uxr)r ] f s e qrs + ^p(Xe lj e pq u!xk)r q fie 3 ki + 

(45) 

£,i£,p{e.ije pq LOxk^Xr)rif s )£jkl£ q rs 

where £ = (^1,^2,^3) an d £ijk = {j — i){k — i){k — j)/2 represents the Levi-Civita tensor arising from 
the cross product. 
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Since A is the average of the velocity increment per unit distance, it is constant with respect the statis- 
tics of Q and u)\, thus (A 2 ...) = A 2 (...). Because of Eq. (|2"2"j) . Q and u>\ are statistically independent, 
so that (e i3 e pq uj\k) = (eije pq ){w\k) =0 and (eije pq uj\kUJ\r} = (e i:j e pq )(uj\ k uj\ r ), therefore second and 
third terms vanish into Eq. (|45[) . 

Due to the isotropy, the Lyapunov basis satisfies Eq. (1231) (X e ij e pql = $ip$jq/3). This is introduced 
in the first term of Eq. (|4"5)) . which thus depends on r 2 alone. This term, which equals A 2 r 2 /3, is 
associated to the direction f (maximal Lyapunov exponent direction) and represents one degree of 
freedom in the space. Again thanks to the isotropy, the last term of Eq. (|45p is two times the first one 
because it is caused by u>\ which corresponds to the two directions orthogonal to f and thus to two 
degrees of freedom. As a result, taking into account that Eijk^ijh = 2$hk, the standard deviation of the 
longitudinal velocity difference is 

{{Au r f) = A 2 r 2 (46) 
with 



AW = V^^ (47) 

This standard deviation can be expressed through the longitudinal correlation function / 

((Au r f) = 2u 2 (l - f(r)) (48) 

being u the standard deviation of the longitudinal velocity. The maximal Lyapunov exponent is calcu- 
lated in function of /, from Eqs. (|4"6")l and (|48l) 



A(r) = -V2(l-/(r)) (49) 

Hence, substituting Eq. (|4Uf into Eq. (|4"4"1) . one obtains the expression of K in terms of / and its 
gradient 



< 50 > 

This is the proposed closure of the von Karman-Howarth equation, the main asset of which on the 
different diffusion models is that it has been derived from a specific Lyapunov analysis, under the 
assumption of reasonable statistical hypotheses about the longitudinal velocity difference. According 
to Eq. (|50[) . if is a nonlinear term of the first order, thus it does not represent a diffusion approximation, 
but rather a nonlinear advection term which makes Eq. (|79[) a nonlinear partial differential equation 
of the first order when v = 0. 

Equation (|50p corresponds to a mechanism of the kinetic energy transfer which preserves the average 
values of the momentum and of the kinetic energy. Specifically, the analytical structure of Ea. (15T)l) 
states that this mechanism consists of a flow of the kinetic energy from large to small scales which 
only redistributes the kinetic energy between wavelengths. 

This mechanism can be interpreted as follows. If, at to, a toroidal volume S(to) is taken which 
contains X and X' (see Fig. [5]), its geometry and position change according to the fluid motion, and 
its dimensions, yfSp and R, vary with the time to preserve the volume. Choosing S in such a way that 
R increases with the time, the Lyapunov analysis of Eqs. (fTBf leads to R rj R(to) e A( - t ~*°' ) . According 
to the theory 20 , for t > to, the trajectories of the two particles are enclosed into S(i). Hence, the 
kinetic energy, initially enclosed into S(to), at the end of the fluctuation is contained into S(t) whose 
dimensions are changed with respect to S(to). The kinetic energy is then transferred far from X and 
X', resulting enclosed in a more thin toroid. 



6 Skewness of the velocity difference PDF 

The obtained expression of K{r) allows to determine the skewness of Au r [5] 

H j r) = <(^> 3 > = m (5i) 

{{Au r ) 2 f 2 (2(l-/(r))) 3/2 
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which is expressed in terms of the longitudinal triple correlation fc(r), linked to K(r) by K(r) = 
u 3 (d/dr + 4/r) k(r) (also see Appendix, Eq. ([511) ). Since / and k are, respectively, even and odd 
functions of r with /(0) = 1, fc(0) = fc'(O) = fc"(0) =0, ff 3 (0) is given by 

^3(0) = limF 3W = I - 77 ^i_ (52 ) 

where the apex denote the derivative with respect to r. To obtain £^(0), observe that, near the origin, 
K behaves as 

K = u 3 ^rW)f"(0)Y + 0(r 4 ) (53) 

then, substituting Eq. ([5"3]) into K(r) — u 3 (d/dr + 4/r) k(r) and accounting for Eq. (j5"2j) . one obtains 

H 3 (0) = ~ = -0.42857... (54) 

H^(0) is a constant of the present analysis, which does not depend on the Reynolds number. This is 
in agreement with the several sources of data existing in the literature such as [5] , [315] , [37] , [35] (and 
Refs. therein) and its value gives the entity of the mechanism of the energy cascade. 

This skewness causes variations in the time of the Taylor scale in accordance to Eq. (JTHJ) ■ The time 
derivative of At is calculated considering the coefficients of the order 0(r 2 ) in the von Karman-Howarth 
equation [7], [8] which are determined substituting Eqs. ((53]) and ([54]) into Eq. (|79l) 



where f IV = d 4 f/dr 4 . The first term tends to decrease At and expresses the mechanism of energy 
cascade, whereas the term in the brackets gives the viscosity effect which, in general, tends to increase 
At depending on the current values of f IV (0) and At [Z], [5]. 



7 Statistical analysis of the velocity difference 

As explained in this section, the Lyapunov analysis of the local deformation and some plausible as- 
sumptions about the statistics of the velocity difference Au(r) = u(X + r) — u(X) lead to determine 
all the statistical moments of Au(r) with only the knowledge of the function K(r) and of the value of 
the critical Reynolds number. 

The statistical properties of Z\u(r), are investigated expressing the velocity fluctuation, given by 
Eq. (fT2"]) . as the Fourier series 

where U(k) = (Ui(k) 1 U2{i^) 1 U 3 (k)) are the components of the velocity spectrum, which satisfy the 
Fourier transformed Navicr-Stokcs equations [8J 

dU p (K) 



Ot 



-vk 2 U p (n) + i YJP^VJ$Ft(k. - j) - K q U q Wri" ~ J)) (57) 



All the components U(k) ~ d\J {n) / dt / A are random variables distributed according to certain distri- 
bution functions, which are statistically orthogonal each other [S]. 

Thanks to the local isotropy, u is sum of several dependent random variables which are identically 
distributed [5J, therefore u tends to a gaussian variable [33J, and U(k) satisfies the Lindeberg condition, 
a very general necessary and sufficient condition for satisfying the central limit theorem [33) . This 
condition does not apply to the Fourier coefficients of Au. In fact, since Au is the difference between 
two dependent gaussian variables, its PDF could be a non gaussian distribution function. In x = 0, 
the velocity difference Z\u(r) = (Aui, Au-2, Au^) is given by 

Au p ^jY,^f^( elKr - 1 ^ L + B + P + N ( 58 ) 
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N = N + Ci - Cl 



This fluctuation consists of the contributions appearing into Eq. (|5T[) : in particular, L represents the 
sum of all linear terms due to the viscosity and B is the sum of all bilinear terms arising from inertia 
and pressure forces. P and N arc, respectively, the sums of definite positive and negative square terms, 
which derive from inertia and pressure forces. The quantity L + B tends to a gaussian random variable 
being the sum of statistically orthogonal terms [39], [33], while P and N do not, as they are linear 
combinations of squares |39) . Their general expressions are |39| 

p = Po + m + nl 

(59) 

■52 

where Po and Nq are constants, and 771, 772, Cl and C2 are four different centered random gaussian 
variables. Therefore, the fluctuation Au r of the longitudinal velocity difference can be written as 

Au p = Vi(r)£ + Mr) (X(V 2 - 1) - (C 2 - 1)) (60) 
where £, rj and C are independent centered random variables which have gaussian distribution functions 
with standard deviation equal to the unity. The parameter x is a positive definite function of the 
Reynolds number, whereas ip% and V>2 are functions of the space coordinates and the Reynolds number. 

At the Kolmogorov scale £, the order of magnitude of the velocity fluctuations is uk 2 t/£, with 
t = 1/A and uk — v/i, whereas is negligible because is due to the inertia forces: this immediately 
identifies ifii « uk 2 t/£. 

On the contrary, at the Taylor scale At, ipi is negligible and the order of magnitude of the velocity 
fluctuations is u 2 t/Xt, therefore 7^2 ~ u 2 t/Xt- 
The ratio ip2 /V'l is a function of R\ 



V'i(r) uk^Xt y 15V15 



where ip(r) = 0(1), is a function which has to be determined. 

Hence, the dimensionless longitudinal velocity difference Au ri is written as 

Au r _ e + ^(x(7? 2 -l)-(C 2 -l)) 



(62) 



The dimensionless statistical moments of Au r are easily calculated considering that £, rj and £ are 
independent gaussian variables 

ff - s *' <f " ><(x( " 2 - 11 - (c2 - m (63) 

where 

k 



{{Au r f) n/2 (l + 2^(i + x 2 )) „/2^ 

((x(v 2 1) - (C 2 - i)) fe > = £ (f) (-x) 4 «C 2 " - !) fe ") 



(64) 



i=0 

In particular, the third moment or skewness, H3, which is responsible for the energy cascade, is 

. ^ (*•-'> (65) 

(1 + 2tA 2 (1 + X 2 )) 3/2 

For x 7^ 1) the skewness and all the odd order moments are different from zero, and for n > 3, all the 
absolute moments are rising functions of R\ , thus Au r exhibits an intermittency whose entity increases 
with the Reynolds number. 

All the statistical moments can be calculated once the function x(R\) and the value of ipo are 
known. The expression of K(r) obtained in the first part of the work allows to identify ^(O) and then 
fixes the relationship between ip and x(-^a) 

-*(»)- ^V 1 -* 3 ' „ -j (60) 

(l + 2*„ 2 (l + ^)) 3 ' 2 1 
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Fig. 3 Parameter \ plotted as the function of R\ 



R, 



where tp = V"(0, R\) = 0{^/R\) and \ = x{R\) > 0. This relationship does not admit solutions with 
X > below a minimum value of (R\) m in dependent on ipg. According to the analysis of section [21 
(R\)min is chosen to 10.12, which corresponds to ipo — 1-075. (setting \ — 0, R\ = 10.12 in ^(O)). 
Varying the value of (Rx)min from 8.5 to 15 would bring values of tpo between 1.2 and 0.9, respectively. 
In figure[31 the function x(R\) is shown for tpo = 1.075. The limit \ ~ 0.86592 for R\ — > oo is reached 
independently of the value of i/'o- 

The PDF of Au r is expressed through the Frobenius-Perron equation 

F(Au'r) = [ [ [ P(0P(V)P(0 S (Au' r - Au r (£, V , Q) d£ dr, d( (67) 
J£ Jn JQ 

where Au r is calculated with Eq. (|62l) . 5 is the Dirac delta and p is a gaussian PDF whose average 
value and standard deviation are equal to and 1, respectively. 

For non-isotropic turbulence or in more complex cases with boundary conditions, the velocity 
spectrum could not satisfy the Lindeberg condition, thus the velocity will be not distrubuted following 
a Gaussian PDF, and Eq. (|60p changes its analytical form and can incorporate more intermittant terms 
33J which give the deviation with respect to the isotropic turbulence. Hence, the absolute statistical 
moments of Au r will be greater than those calculated with Eq. (|52l . indicating that, in a more complex 
situation than the isotropic turbulence, the intermittency of Au r can be significantly stronger. 



8 Results and discussion 

In order to validate the results of the proposed Lyapunov analysis, several data are now presented. 

As the first result, the evolution in the time of / is calculated with the proposed closure (Eq. 
([50)1 ) of the von Karman-Howarth equation, where the boundary conditions are given by Eq. ([52"]) (see 
Appendix). The turbulent kinetic energy and E{n) and T(k) are calculated with Eq. ([55)1 and Eqs. 
([Mf . respectively. The calculation is carried out for an initial Reynolds number of Re(0) — u{i))L r /v 
= 2000, where L r and u(0) are, respectively, the reference dimension and the initial velocity standard 
deviation. The initial condition for / is f(r) = cxp (— l/2(r/Ar) 2 ), where Ar/L r — l/(2-\/2), whereas 
u(0) — 1. The dimensionless time of the problem is defined as i — t u(0)/L r . 

Equation (f79"|) was numerically solved adopting the Crank- Nicholson integrator scheme with variable 
time step At, whereas the discretization in the spatial domain is made by N — 1 intervals of the same 
amplitude Ar. The truncation error of the scheme of integration is 0(At 2 ) + 0(Ar 2 ), where At is 
automatically selected by the algorithm in such a way that 0(At 2 ) <~ 0(Ar 2 ) for each time step, 
thus the accuracy of the scheme is of the order of 0(Ar 2 ). As the consequence, the discretization of 
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the Fourier space is made by N — 1 subsets in the interval [0,km], where km — ir/(2Ar). For the 
adopted initial Reynolds number the choice N — 1500 gives an adequate discretization, which provides 
Ar < I, for the whole simulation. For what concerns u, it was calculated with Eq. (I55|) and the kinetic 
energy was checked to be equal to the integral over k of E(k). During the simulation, T(k) must 
identically satisfy Eg. (1551) (see Appendix) which states that T(k) does not modify the kinetic energy. 
The integral of T(k) is calculated with the trapezoidal formula, for k £ (0, km), arid the simulation 
will be considered to be accurate as long as 



T(k)cIk - 



T(k)cLk = 



(68) 



namely, when T{k) ~ for k > km- As the simulation advances, according to Eq. (1501) . the energy 
cascade determines variations of E{k) and T(k) for values of k which rise with the time and that can 
occur out of the interval (0, km)- Thus, Eq. Q68[) holds until a certain time, where these wave-numbers 
are about equal to km- For higher times, the variations of T(k) occur for k > km, out of the interval of 
numerical integration (0, km), and Eq. (|68p could not be satisfied. Hence, the results of the simulation 
will be accurate until reaching of the following condition [4"0] 



T(k)cIk 



> 



1 K 



M 



12 N 2 



d 2 T 



Ok 2 



(69) 



where the right-hand side represents the estimation of the truncation error of the trapezoid formula 
40 . There, it is found that, Ar w 0.8 £ in all the simulations. . Thereafter, the numerical method loses 
its consistency, since the violation of Eq. (|6"8")) implies that K does not preserve the kinetic energy. 
Nevertheless, the error committed by the algorithm is considered still acceptable as long as [7], [5] 



< 



T{k)(Lk 



dt 



(70) 
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Fig. 5 Plot of E(k) and T(k), for Re(0) = 2000, at the diverse times of simulation. 

That is to say, when the numerical residual of the rate of energy transfer is much less than the 
dissipation rate. In any case the simulations are stopped if £ reaches its relative minimum in function 
of the time with £ > Ar, or as soon as Ar = £. The accuracy of the solutions is evaluated through the 
check of Eq. (|70p . At end of simulation, we assume that the energy spectrum is fully developed. 

In order to study the proposed closure, we first analyze the case with K = 0. In this case, with the 
assumed initial condition, Eq. (|79l) admits the analytical self-similar solution [7] 

f(t, r) = cxp (~ (x^tj) ) ' Wh6re AtW = \At(°)+ 4 ^ ( 71 ) 

f(t,r) maintains its shape unchanged in function of the dimensionless coordinate r/Ar(i)- The corre- 
sponding numerical solutions were calculated for different spatial discretization (i.e. N = 500, 1000, 
and 1500). All these simulations satisfy Eq. (|7ip . with a calculated error which is always of the order 
of the truncation error of the scheme of integration. Also the energy spectrum preserves its shape, and 
is given by E(k\ t ) = Au 2 (kX t ) A exp(-a(KA T ) 2 ), with A = 0(1), a = 0(1) [8]. 

Consider now the case, where K is given by Eq. (|5U)) . The diagrams of Fig. [U show the evolution 
of f(r) and k(r) in terms of r/Xx, at different times of simulation. The kinetic energy and At vary 
according to Eqs. (|50l) . ([83]) and (j55|) . thus f(r) and k(r) change in such a way that their scales, in 
particular £ and At, diminish as the time increases, whereas the maximum of |fc| decreases. That is, 
Eq. (|50[) corresponds to a transferring of the energy toward the smaller scales which strongly contrasts 
the effects of viscosity seen in the previous case. 

At the final instants of the simulation, one obtains that / — 1 = 0( r 2 / 3 ) for r/Xr — 0(1), and the 
maximum of |fc| is about 0.05. These results are in very good agreement with the numerous data of the 
literature ;8j which concern the evolution of / in homogeneous isotropic turbulence. Figure [S] shows 
the diagrams of E(n) and T(k), at the same times, in comparison with the Kolmogorov law (k~ 5 / 3 ) 
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Fig. 6 Plot of the final energy spectrum for different spatial discretization, Re(0) = 2000. 

and with the incompressibility condition (k 4 ). The energy spectrum depends on the initial condition, 
and at the end of the simulation, can be compared with the Kolmogorov spectrum in an opportune 
interval of wave-numbers which defines the inertial subrange. This arises from /, which, at the final 
times, behaves like / — 1 = O (r 2 / 3 ) for r = O(Xt)- E{k) satisfies the continuity equation which 
implies that E{k) sa k 4 near the origin. In the figure, the dimensionless time t = 0.63 corresponds to 
the condition whereas i = 0.69 (bold dashdotted line) represents the spectrum calculated when 
£ reaches its minimum in function of the time. In this last situation, the algorithm is less accurate, 
due to the discretization. In fact, T(k), (bold dashdotted), exhibits small nonzero values in a wide 
range of variations of k, for k > km- Although Eq. (1681) is not satisfied, the error of the algorithm 
is still acceptable since the absolute value of the integral of T(«) over (0, km) is about four orders of 
magnitude less than the dissipation rate (see Eq. ((70)) '). 

To study the effect of the spatial discretization on the final energy spectrum, Fig. [5] shows the 
results of other simulations which were performed for N — 500, 1000. It is found that, for N = 500 
and 1000, the stopping criterion is satisfied for I = Ar, whereas the final times of simulation rises with 
N, resulting that i w 0.33 and i ~ 0.64 for N = 500 and 1000, respectively. As the consequence, also 
the inertial subrange of Kolmogorov increase with N. 

In order to compare the results of the present analysis with the data in the literature, a further 
calculation has been carried out with an initial Reynolds number of Re(0) — 3000. The present results 
(continuous line) are shown in Fig. [3 in terms of energy spectrum, and are compared with those 
calculated, for the same conditions, with the Oberlack's model [12] (dashed line). For this latter, the 
parameter fc 2 HU (see introduction) is assumed to be equal to k 2 = ~H 3 (0)/(^/2 6). For this value of 
k 2 , the Oberlack's model gives the same value of skewness H 3 (0) = —3/7, calculated with the present 
analysis. Thus, the entity of the mechanism of the energy cascade is the same in both the cases. Again, 
the initial correlation function is / = exp(— (r/Ar) 2 /2), Xt / L r — l/(2\/2) and N = 1500. 

Equation (l68l) is satisfied until a certain time (i.e. t ~ 0.66 for Oberlack's model and t 0.62 for 
the present analysis), therafter the numerical scheme exhibits a minor accuracy, since T(k) ^ for 
k > km)- In the figure, the energy spectrum is shown at the end of the two simulations which happen at 
t = 0.71 (Oberlack) and at i = 0.67 (Present Analysis), where the stopping condition is satisfied for £ = 
min rj Ar. There, the integral of T(k) is, in both the cases, much less than \du 2 /dt\ (about four orders 
of magnitude) and for this reason the results of the two simulations can be considered accurate enough. 
Since the skewness H^{Q) is the same in the two cases, the variations of At almost coincide during 
the simulations, a part very small variations caused by v (see Eq. (I55[) ) which are not appreciated 
in the figure. Therefore, also du 2 /dt is about the same in both the cases. Nevertheless, the energy 
spectra show significative differences. The spectrum calculated with Eq. (I50p exhibits a wider range 
of wave-numbers than the other one, and this is due to the fact that Eq. ([5U)) represents a nonlinear 
term of the first order which does not cause any diffusion effect, whereas for the Oberlack's model, the 
closure term is of the second order and the variable eddy diffusivity produces a sizable reduction of all 
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Fig. 7 Comparison of the results for Re(0)=3000: Continuous line: present analysis. Dashed line: Oberlack's 
model (Ref. 12.). (a) Dotted curve: initial condition, (b) Kinetic energy, (c) Taylor scale 



the wave-numbers, especially for large k. The wave-number calculated where E(k) — 10~ 6 (see figure), 
is quite different in the two cases. Specifically, the present analysis gives a value two times greater than 
in the Oberlack's model, and, as the consequence, also the Kolmogorov inertial subrange calculated 
with the present analysis is about two times the other one. 

Next, the Kolmogorov function Q(r) and Kolmogorov constant C, are determined, using the results 
of the previous simulations. 

Following the Kolmogorov theory, the Kolmogorov function, which is defined as 



Q(f) 



re 



(72) 



is constant w. r. t. r, and is equal to 4/5 as long as t/Xt = 0(1). As shown in Fig. [51 for t — 0, Q m ax 
is significantly greater than 4/5 and the variations of Q with t/Xt can not be neglected. This is the 
consequence of the choice of the initial correlation function. At the successive times, Q ma x decreases 
until the final instants, where, with the exception of r /At ~ 0, Q{r) exhibits variations which are less 
than those calculated at the previous times in a wide range of r/Ay, with a maximum which can be 
compared to 0.8. The Kolmogorov constant C is also calculated by 



C= max [^-) (73) 

For Re(0) = 2000 and t ~ 0.63, C ~ 1.932and Q max ~ .73, whereas for t ~ 0.69, C ~ 1.92, Q max ~ .72. 
For Re(0) — 3000 at end simulation (i.e. t ~ 0.67), C ~ 1.941, Q max — -75, namely C and Q max agree 
with the corresponding quantities known from the literature. 

For Re(0) = 2000, Fig. shows the maximal finite scale Lyapunov exponent, calculated with Eq. 
(|4"9"| . For t = 0, the variations of A are the result of the adopted initial correlation function which is a 



5/3 
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gaussian, whereas as t increases, the variations of / determine sizable increments of A and of its slope 
in proximity of the origin. Then, for t — 0.6, since / — 1 ~ 0(r 2 / 3 ), the maximal finite scale Lyapunov 
exponent behaves like A ~ r~ 2 / 3 . Thus, the Richardson's diffusivity associated to the relative motion 
between two fluid particles, defined as Dr oc Xr 2 [5], here satisfies the famous Richardson scaling law 
Dr sa r 4 / 3 [S] in proximity of the end of the simulation. 

In the figures lit) and [St, skewness and flatness of Au r are shown in terms of r for t — and 0.6 
and Re(0) = 2000. The skewness, H3 is first calculated with Eq. (f5T|) . then H4 has been determined 
using Eq. (j6"3")l . At i = 0, \Ha\ starts from 3/7 at the origin with small slope, then decreases until 
reaching small values. H4 also exhibits small derivatives near the origin, where H4 3> 3, thereafter it 
decreases more rapidly than |i?3|. At t =0.6, the diagram importantly changes and exhibits different 
shapes. The Taylor scale and R\ are both changed, so that the variations of H3 and H4 are associated 
to smaller distances, whereas the flatness at the origin is slightly less than that at t — 0. Nevertheless, 
these variations correspond to higher values of r/Xr than those for t = 0, and again H4 reaches the 
value of 3 more rapidly than H3 tends to zero. 

The PDFs of Au r are calculated, for Re(0) = 2000, with Eqs. and (|52"|). and are shown in Fig. 
1101 in terms of the dimensionless abscissa 

Au r 
S = ((Zki r ) 2 ) 1/2 

where, these distribution functions are normalized, in order that their standard deviations are equal 
to the unity. The figure represents the distribution functions of s for several t/Xt, at i — 0, 0.5 and 
0.6, where the dotted curves represent the gaussian distribution functions. The calculation of H^{r) 
is first carried out with Eq. (|5"Tj) . then ip(r,R\) is identified through Eq. (|6"5"1) . and finally the PDF is 
obtained with Eq. (|57]). For t = (see Fig. [TUh ) and according to the evolutions of H 3 and i?4, the 
PDFs calculated at r/Xr — and 1, are quite similar each other, whereas for r/Xx = 5, the PDF 
is almost a gaussian function. Toward the end of the simulation, (see Fig. [TUb and c) , the two PDFs 
calculated at r/Xr = and 1, exhibit more sizable differences, whereas for r/Xr — 5, the PDF differs 
very much from a gaussian PDF. This is in line with the plots of H%{r) and H^ir) of Fig. |H1 

Next, the spatial structure of Au r , given by Eq. (|B"2"|) . is analyzed using the previous results. Ac- 
cording to the various works |41j . [12"] . Au r behaves quite similarly to a multifractal system, where 
Au r obeys to a law of the kind Au r (r) r q where the exponent q is a fluctuating function of the 
space coordinates. This implies that the statistical moments of Au r (r) are expressed through different 
scaling exponents C(P) whose values depend on the moment order P, i.e. 

{(Au r ) p (r)) - A P r<( p ) (74) 
These scaling exponents are here identified through a best fitting procedure, in the interval (a, a 
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+At)j where the endpoints a is an unknown quantity which has to be determined. The location of 
this interval varies with the time. The calculation of a, £p and Ap is carried out through a minimum 
square method which, for each moment order, is applied to the following optimization problem 

M(P,A P )= (((Au r f) - Apr^fdr = min, P= 1,2,... (75) 

J a 

where ((Au^)) are calculated with Eqs. (j53"]l . and a is calculated in order to obtain £(3) = 1. 

Figure [11] shows the comparison between the scaling exponents here obtained (continuous lines with 

solid symbols) and those of the Kolmogorov theories K41 (dashed lines) and K62 [3T] (dashdotted 

Table 1 Scaling exponents of the longitudinal velocity difference for different initial Reynolds numbers. 
Re(0) P 1 2 3 4 5 6 7 8 9 10 11 12 13 14 



2000 C(P) 0.36 0.71 1.00 1.19 1.41 1.61 1.84 2.04 2.25 2.49 2.72 2.93 3.15 3.38 
3000 C(P) 0-36 0.72 1.00 1.19 1.43 1.64 1.84 2.03 2.25 2.49 2.71 2.92 3.11 3.33 
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Fig. 10 PDF of the velocity difference fluctuations at the times t —0 (a), t = 0.5 (b) and t —0.6 (c). Continuous 
lines are for r —0, dashed lines are for r/Xr = 1, dot-dashed lines are for r/Xr =5, dotted lines are for gaussian 
PDF. 



lines), and those given by She-Leveque [H] (dotted curves). At t = 0, the values of C(P) are the result 
of the chosen initial condition. As the time increases, / changes causing variations in the statistical 
moments of Au r (r). As result, £(P) gradually diminish and exhibit a variable slope which depends on 
the moment order P, until to reach the situation of Fig. [TTb . where the simulation is just ended. The 
dimensionless moments of Au r (r) are changed. The plot of C(P) shows that near the origin, ((P) ~ -P/3, 
and that the values of C(P) seem to be in agreement with the those proposed by She-Leveque. More in 
detail, Table I reports these scaling exponents in terms of the moments order, calculated for t = 0.63. 
These values and the peculiar diagram of Fig. (flT]) are the consequence of the spatial variations of 
the skewness, calculated using Eq. (|5ip . and of the quadratic terms due to the inertia and pressure 
forces into the expression of the velocity difference, which make ((Au r ) p ) a quantity quite similar to 
a multifractal system. 

Other simulations with different initial correlation functions and Reynolds numbers have been 
carried out, and all of them lead to analogous results, in the sense that, at the end of the simulations, 
the diverse quantities such as Q(r), C and C(P) are quite similar to those just calculated. For what 
concerns the effect of the Reynolds number, its increment determines a wider range of the wave- 
numbers where E(k) is comparable with the Kolmogorov law and a smaller dissipation energy rate in 
accordance to Eq. 

In order to study the evolution of the intermittency vs. the Reynolds number, Table II gives the 
first ten statistical moments of F(du r /dr). These are calculated with Eqs. (|63|) and (j64|) . for R\ = 
10.12, 100 and 1000, and are shown in comparison with those of a gaussian distribution function. It is 
apparent that a constant nonzero skewness of du r /dr, causes an intermittency which rises with R\ (see 
Eq. (jni])). More specifically, Fig. [T2l shows the variations of #4(0) and Hq(0) (continuous lines) in terms 
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t=0.63 




Fig. 11 Scaling exponents of the longitudinal velocity difference versus the order moment at different times, 
Re(0) — 2000. Continuous lines with solid symbols are for th e pr esent data. Dashed lines are for Kolmogorov 
K41 data [6]. Dashdotted lines are for Kolmogorov K62 data [41] . Dotted lines are for She-Leveque data [42] 



of R\, calculated with Eqs. (|63p and with i?3(0) = —3/7. These moments are rising functions 

of R\ for 10 < R\ < 700, whereas for higher R\ these tend to the saturation and such behavior 
also happens for the other absolute moments. According to Eq. (jf>3"]) . in the interval 10 < R\ < 70, 
H4 and Hq result to be about proportional to _R° 34 and -R°' 78 , respectively, and the intermittency 
increases with the Reynolds number until R\ ~ 700, where it ceases to rise so quickly. This behavior, 
represented by the continuous lines, depends on the fact that ip ~ V R\, and results to be in very good 
agreement with the data of Pullin and Saffman [33]; for 10 < R\ < 100. Figure can be compared 
with the data collected by Sreenivasan and Antonia [35], which are here reported into Fig. [T3J These 
latter are referred to several measurements and simulations obtained in different situations which can 



Table 2 Dimensionless statistical moments of F(du r /dr) at different Taylor scale Reynolds numbers. P.R. as 
for "present results". 

Moment R x « 10 R x = W R x = 10 3 Gaussian 
Order P. R. P. R. P. R. Moment 



3 


-.428571 


-.428571 


-.428571 





4 


3.96973 


7.69530 


8.95525 


3 


5 


-7.21043 


-11.7922 


-12.7656 





6 


42.4092 


173.992 


228.486 


15 


7 


-170.850 


-551.972 


-667.237 





8 


1035.22 


7968.33 


11648.2 


105 


9 


-6329.64 


-41477.9 


-56151.4 





10 


45632.5 


617583. 


997938. 


945 
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Fig. 12 Dimensionless moments i?4(0) and He(0) plotted vs. J?a- Continuous lines are for the present results. 
The dashed line is the tangent to the curve of -ff.i(0) in R\ ~ 10. 
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Fig. 13 Flatness i/ 4 (0) vs. R\. These data are from Ref. [38] 



be very far from the isotropy and homogeneity conditions. Nevertheless a comparison between the 
present results and those of Ref. [35] is an opportunity to state if the two data exhibit elements in 
common. According to Ref. [38] . the flatness monotonically rises with R\ with a rising rate which 
agrees with Eq. (|64l) for 10 < R\ < 60 (dashed line, Fig. 1X2]) . whereas the skewness seems to exhibit 
minor variations. Thereafter, Hi continues to rise with about the same rate, without the saturation 
observed in Fig. [T^l The weaker intermittency calculated with the present analysis seems arise from 
the isotropy which makes u r a gaussian random variable, while, as seen in sec. [SI without the isotropy, 
the flatness of u r and Au r can be much greater than that of the isotropic case. 

Next, the obtained results are compared with the data of Tabeling et al [3J>], [SZ]- There, in an 
experiment using low temperature helium gas between two counter-rotating cylinders (closed cell), 
the authors measure the PDF of du r /dr and its moments. Again, the flow could be quite far from 
to the isotropy condition, since these experiments pertain wall-bounded flows, where the walls could 
importantly influence the fluid velocity in proximity of the probe. The authors found that moments 
Hp, with p > 3, first increase with R\ until R\ ss 700, then exhibit a lightly non-monotonic evolution 
with respect to R\, and finally cease their variations, denoting a transition behavior (See Fig. 
As far as the skewness is concerned, the authors observe small percentage variations. Although the 
isotropy does not describe the non-monotonic evolution near R\ = 700, the results obtained with Eq. 
(|62| can be considered comparable with those of Refs. [55], [37], resulting again, that the proposed 
analysis gives a weaker intermittency with respect to Refs. [33], |37j . 
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s 

Fig. 15 Log linear plot of the PDF of du r /dr for different Rx- (a): dotted, dashdotted and continuous lines 
are for R\ = 15, 30 and 60, respectively, (b) and (c) PDFs for R\ = 255, 416, 514, 1035 and 1553. (c) represents 
an enlarged part of the diagram (b) 



The normalized PDFs of du r /dr are calculated with Eqs. (|6"T|) and (|62p . and are shown in Fig. (TS] 
in terms of the variable s, which is defined as 

du r /dr 

*~ {{dur/dr) 2 ) 1 ' 2 

Figure IT5a shows the diagrams for R\ — 15, 30 and 60, where the PDFs vary in such a way that 
# 3 (0) = -3/7. 

As well as in Ref. [37J, Figs. 4b and 4c give the PDF for R x = 255, 416, 514, 1035 and 1553, where 
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Fig. 16 PDF of dur/dr for R x = 255, 416, 514, 1035 and 1553. These data are from Ref. [37] 
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Fig. 17 Plot of the integrand s 4 F(s) for different R\. Dotted, dashdotted and continuous lines are for R\ = 
15, 30 and 60, respectively. 

these last Reynolds numbers are calculated through the Kolmogorov function given in Ref. [37] , with 
-^3(0) = —3/7. In particular, Fig. [T5b represents the enlarged region of Fig. [T5b . where the tails of 
PDF are shown for 5 < s < 8. According to Eq. the tails of the PDF rise in the interval 10 

< R\ < 700, whereas at higher values of R\, smaller variations occur. Although the trend observed in 
Fig.Ql)] 37] is non-monotonic, Fig.lPob shows that the values of the PDFs calculated with the proposed 
analysis, for 5 < s < 8, exhibit the same order of magnitude of those obtained by Tabeling et al [37]. 

Asymmetry and intermittency of the distribution functions are also represented through the inte- 
grand function of the 4 th order moment of PDF, which is J^s) = s 4 F(s). This function is shown in 
terms of s, in Fig. (T7J for R\ = 15, 30 and 60. 



9 Conclusions 

The proposed analysis is based on the conjecture that the turbulence is caused by the bifurcations 
and that the kinematic of the relative motion is much more rapid than the fluid state variables. 
The analysis also assumes statistical hypotheses about the velocity difference, which derive from the 
condition of fully developed turbulence. The main limitation of this analysis is that it only studies the 
developed homogeneous-isotropic turbulence, whereas this does not consider the intermediate stages 
of the turbulence. 

The results can be here summarized: 
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1. The qualitative analysis of the bifurcations leads to determine the order of magnitude of the critical 
Reynolds number based on the Taylor scale and the number of the bifurcations at the onset of the 
turbulence. 

2. The momentum equations written using the referential coordinates allow to factorize the velocity 
fluctuation and to express it in Lyapunov exponential form of the local fluid deformation. As a 
result, the velocity fluctuation is the combined effect of the exponential growth rate and of the 
rotations of the Lyapunov basis with respect to the fixed frame of reference. 

3. The finite scale Lyapunov analysis of the relative motion provides an explanation of the physical 
mechanism of the energy cascade in turbulence and gives a closure of the von Karman-Howarth 
equation. This is a non-diffusive closure which depends on the local values of the longitudinal 
correlation function and on its spatial gradient. 

The fluid incompressibility is a sufficient condition to state that the inertia forces transfer the 
kinetic energy between the length scales without changing the total kinetic energy. This implies 
that the skewness of the longitudinal velocity derivative is a constant of the present analysis and 
that the energy cascade mechanism does not depend on the Reynolds number. 

4. The statistics of Au r can be inferred looking at the Fourier series of the velocity difference. This 
is a non-Gaussian statistics, where the constant skewness of du r /dr implies that the other higher 
absolute moments increase with the Taylor-scale Reynolds number. 

5. The proposed closure of the von Karman-Howarth equation, shows that the mechanism of en- 
ergy cascade gives energy spectrum that can be compared with the Kolmogorov law k -5 / 3 in an 
opportune range of wave-numbers and which satisfy the incompressibility condition. 

6. For developed energy spectrum, the Kolmogorov function exhibits, in an opportune range of r, 
small variations much less than at the previous times, and its maximum is quite close to 4/5, 
whereas the Kolmogorov constant is about equal to 1.93. As the consequence, the maximal finite 
scale Lyapunov exponent and the diffusivity coefficient vary according to the Richardson law when 
the separation distance is of the order of the Taylor scale. 

7. The analysis also determines the scaling exponents of the moments of the longitudinal velocity 
difference through a best fitting procedure. For developed energy spectrum, these exponents show 
variations with the moment order consistent with those present in the literature. 



10 Appendix 

The von Karman-Howarth equation gives the evolution in the time of the longitudinal correlation 
function for isotropic turbulence. The correlation function of the velocity components is the symmetrical 
second order tensor Rij(r) = (uiUj}, where Ui and u'j are the velocity components at x and x + r, 
respectively, being r the separation vector. The equations for i£y are obtained by the Navier-Stokes 
equations written in the two points x and x + r [7] , [5] . For isotropic turbulence Rij can be expressed 
as 

R tJ (r) = u 2 [(f-g)^+gS ll ] (76) 

/ and g are, respectively, longitudinal and lateral correlation functions, which are 

(u r (xK(x + r)) (u„(r)u ra (x + r)) 
f(r) = 2 , g(r) = - % (77 

where u r and u n are, respectively, the velocity components parallel and normal to r, whereas r = |r| 
and u 2 = (u 2 r ) ={ul)= 1/3 ( UiUi). Due to the continuity equation, / and g are linked each other by 
the relationship 

9 = f + \%r (78) 

The von Karman-Howarth equation reads as follows [7J , [5] 

df_K 2v fd 2 f 4<9/\ 1Q „ d2 f (Q)f (7Q) 
dt u 2 \ dr 2 r dr J dr 2 
where K is an even function of r, defined as [7J, [5] 

(r^ + 3] K{r) = A (uXiu, - u' k )) (80) 
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and which can also be expressed in terms of the longitudinal triple correlation function k(r) 
(i^(x)w r (x + r))/ii 3 



K(r) -u^l + J) k(r) 

The boundary conditions of Eq. ([79)) are [7] , [8] 
/(0) = 1, lim f(r) = 



(81) 



(82) 



(83) 



The viscosity is responsible for the decay of the turbulent kinetic energy, the rate of which is [7J, [5] 

This energy is distributed at different wave-lengths according to the energy spectrum E(k) which 
is calculated as the Fourier Transform of fu 2 , whereas the "transfer function" T(k) is the Fourier 
Transform of K [8], i.e. 



'E{k)~ 


i r°° 


~u 2 f(r)~ 


T{k)_ 




K{r) 



2 2 

k r 



cos kt dr 



(84) 



(85) 



where k, = |rc| and T(k) identically satisfies to the integral condition 

I T(n)dn = 
Jo 

which states that K does not modify the total kinetic energy. The rate of energy dissipation e is 
calculated for isotropic turbulence as follows [5] 

,2 foo 

(86) 



3d" 2 „ f °° 2r , v , 

— = 2v I K 2 EU)dn 

2 dt J V ; 



The microscales of Taylor At, and of Kolmogorov £, are defined as 



1 



((dur/dr) 2 ) 9 2 //ar 2 (0)' 



3\ V4 



(87) 
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